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1. Motivation and objectives 

The Smagorinsky subgrid-scale eddy viscosity model (SGS-EVM) is commonly 
used in large eddy simulations (LES) to represent the effects of the unresolved scales 
on the resolved scales. This model is known to be limited because its constant must 
be optimized in different flows, and it must be modified with a damping function to 
, account for near- wall effects. The recent dynamic model (Germano et al. 1991) is 
designed to overcome these limitations but is compositionally intensive as compared 
to the traditional SGS-EVM. In a recent study using direct numerical simulation 
data, Horiuti (1993) has shown that these drawbacks are due mainly to the use 
of an improper velocity scale in the SGS-EVM. He also proposed the use of the 
subgrid-scale normal stress as a new velocity scale that was inspired by a high-order 
anisotropic representation model (Horiuti 1990). The testing of Horiuti (1993), 
however, was conducted using DNS data from a low Reynolds number channel flow 
’ simulation. It was felt that further testing at higher Reynolds numbers and also 
using different flows (other than wall-bounded shear flows) were necessary steps 
needed to establish the validity of the new model. This is the primary motivation 
of the present study. The objective is to test the new model using DNS databases 
of high Reynolds number channel and fully developed turbulent mixing layer flows. 
The use of both channel (wall- bounded) and mixing layer flows is important for 
the development of accurate LES models because these two flows encompass many 
characteristic features of complex turbulent flows. 

2. Accomplishments 

The subgrid-scale stress tensor, T{j , that results from filtering the Navier-Stokes 
equations consists of three terms (Bardina 1983): 

Tij = Lij + Cij + Rij, (1) 


L t j = U t Uj — UiUj, Cij = Itillj + l/'-Uj, Rij = u[uj 

where overlineui denotes the filtered velocity component and u\ = Ui — Uj denotes 
the SGS component of u,. L t j is the Leonard term, C tJ is the cross term, and R tJ 
is the SGS Reynolds stress. The indices i = 1,2,3 correspond to the directions 
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x, y, and z, respectively, with x the streamwise (it, — u), y the “major-gradient 
(wall-normal or cross-stream) (u 2 = v), and 2 the spanwise (u 3 = w) directions. 

The Leonard term in eq. (1) is not modeled but is treated explicitly by applying 
the filter, while the other two terms (C,y and R, } ) need to be modeled. A successful 
model for the cross term is a model suggested by Bardina (1983) where 


Cij = u'iUj + UiU'j 


This model has been tested by Bardina (1983) for homogeneous flows and by Horiuti 
(1989) for the channel flow and was found to be a good model for the cross terms. 
This model will not be tested further in this work. 

For the Rij terms, the eddy viscosity model by Smagorinsky (Smagorinsky 1963): 


2— f ,dtii dxij 


( 2 ) 


I / e = (C 5 A) 2 (i5 0 5, i ] 1 / 2 , 


_ du, dUj 
Sij 


and the Bardina model _ _ 

Rij ~ C(ui - - Uj). (3) 

are two of several models which are used in L ES computations. In these models, 
C s and C are model constants, and E G = u\u\/2 and v e are, respectively, the SGS 
turbulent kinetic energy and SGS eddy viscosity coefficient. A is the characteristic 
SGS length scale whose value is defined as (Aa*AyAz) 1/3 ; Ax, Ay, and Az are the 
grid intervals in the x, y, and z directions, respectively. The Smagorinski model 
is a “Prandtl-type” mixing length model that can be derived by starting with the 
eddy viscosity approximation to the subgrid-scale Reynolds stresses and assuming 
production and dissipation are in balance. In an eddy viscosity approximation, v t 
is written as the product of a characteristic time scale r and a velocity scale 


v t = C u rE 


(4) 


where C v is a model constant, r is then expressed as (Horiuti 1993) 


r = 



dv\dv\ 
dx\ dx\ 



(5) 


where e is the dissipation rate of E G &nd C € is a model constant. The Smagorinsky 
model assumes that E ~ Eq in (4). 

In the present study, we make use of the direct numerical simulation flow fields 
available at CTR to directly test the various approximations. The fields we con- 
sider are homogeneous in two- directions. To compute the large-eddy flow fields, we 
filter the DNS fields by applying a two-dimensional Gaussian filter in the i — 1,3 
directions. In the inhomogeneous direction ( i = 2), a top-hat filter is applied to the 
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FIGURE 1. y-distribution of the SGS-Reynolds shear stress in turbulent channel at 
Re r = 790. (Model with E = E G /u 2 T .) 



FIGURE 2. y-distribution of the SGS-Reynolds shear stress in turbulent channel at 
Re r = 790. (Model with E = u^u^/u 2 .) 

channel flow fields. No filter was applied in this direction ( i = 2) to the mixing layer 
flow field. This is due to the fact that occasionally the doubly filtered ([]) grid-scale 
variables were larger than the singly filtered ones ([]), owing to the inaccuracy of a 
top-hat filter in regions where grid spacing is coarse. 

The DNS databases we used were the fully developed incompressible channel flows 
at Re T (Reynolds number based on the wall-friction velocity, u r , and the channel 
height, 26)= 360 (Kim et al. 1987) and 790 (Kim 1990), and the incompressible 
mixing layer at Reg (the Reynolds number based on the momentum thickness, 6 m , 
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and the velocity difference, A U)= 2400 (Moser and Rogers 1992). We started with 
the low Reynolds number channel flow data as a confidence test. We found that 
the results obtained in this case are consistent with the previous work of Horiuti 
(1993), who used a different set of DNS data but at the same Reynolds number. 
The details of this testing are not shown in the present report. 

The high Reynolds number channel flow field (with 256 x 193 x 192 grid points) was 
filtered to 64 x 97 x 48 grid points. The mixing layer flow field (with 512 x 210 x 192 
grid points) was filtered to 64 x 210 x 48 grid points. These LES grid point numbers 
were chosen so that the turbulent kinetic energy retained in the SGS components 
is large. This is needed to make a fair assessment of the SGS models. SGS model 
evaluations were conducted by comparing the y-distribution of the mean values 
averaged in the x — z plane (denoted by <()>), and also by comparing the 
y-distribution of the root-mean-square (rms) values of the exact terms with the 
model predictions. Only the y-distribution of the mean values are shown in the 
present report because the rms values were found to give similar results. 

2.1 A proper eddy viscosity velocity scale 


2.1.1 Channel flow 

The y-distribution of the SGS Reynolds shear stress < u\u 2 > obtained with 
E = Eg and C u = 0.1 in (4) is compared with the DNS data in Fig. 1. While 
the agreement between the model and the term is good in the central portion of 
the channel, the agreement deteriorates near the wall where the model predicts 
a very large peak compared to the actual data. This overprediction of the shear 
stress near the wall when Eg is used for E in (4) implies that a damping function 
is needed to account for the presence of the wall. This near-wall overprediction of 
the stress is similar to the near- wall behavior of one-point closure models (see Rodi 
& Mansour 1991). This behavior of one-point closure models is attributed to the 
rapid reduction of the Reynolds shear stress (as the wall is approached) due to the 
preferential damping of the normal stress (Launder 1987, and Durbin 1992). Horiuti 
(1993) reasoned that the same wall dampi ng ef fects should hold true for the SGS 
field. Indeed, when the SGS normal stress u 2 u 2 is used for E (with C v = 0.23, see 
Fig. 2), the model agrees well with SGS Reynolds shear stress near the wall without 
an additional damping function. The model is, however, less effective a s compared 
to using the total energy in the core region of the channel. The main deficiency in 
the core region is attributed to excessive grid stretching in the y-direction because 
of the mapping used in conjunction with Chebyshev expansions. In an actual LES 
computation, finite differences with a more uniform grid are used in the y direction 
and, therefore, a more isotropic energy distribution can be expected in this case. 

The effects of the anisotropic grid can be evidenced by the y-distribution profile 
of the ‘flatness parameter’ A (Lumley 1978) averaged in the x-z plane. In this case 
A is defined as 


A = [1 - -{^2 - A$ } 1 , .4 2 = ctijCLij, A 3 = a t jajkaki, 
o 


( 6 ) 
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FIGURE 3. y-distribution of the flatness parameter A and the Van Driest function 
(channel at Re r = 790). 

dij = -^iju' k u' k } /k,k = i u' k u' k 

We find (see Fig. 3) that in the core region of the channel, A ~ 0.35, which is much 
smaller than the expected A = 1 when the small scale turbulence is isotropic. In 
the region around y ~ 0.1, A peaks around A ~ 0.5, and then gradually decreases 
to 0.35 at the channel center. The y-distribution of A for the unfiltered DNS data 
does not show this overshoot and is close to A = 1 around the centerline. The 
grid spacing in the central region of channel seems to be too coarse; therefore, a 
considerable anisotropy exists in the S GS turbulence fluctuations. In fact, when 
the SGS-EVM model with E = u' 2 u\ in (4) was used in an actual LES channel 
flow calculations using a more uniform grid at high Reynolds number (Re r = 1280) 
(Horiuti 1993), a good agreement with experimental data was found. The present 
comparisons for the high Re channel flow confirm the conclusions of Horiuti (1993) 
based on the low Re channel flow fields. 

For the record, the y-distribution of the conventional Van Driest damping func- 
tion ((1 — exp(— y + /26.0)) (normalized with value of A at the channel center) is 
included in Fig. 3. It should be noted that the ‘flatness parameter’ A has a similar 
distribution across the channel as the Van Driest function, suggesting that A may 
be used as an alternative method to damp the eddy viscosity near the wall (Horiuti 
1992). 

2.1.2 Mixing layer 

The y-distribution of u[u 2 obtained using E = Eg ( C u = 0.20) and E = u' 2 u 2 
( C v = 0.26) in (4) are compared with the DNS data in Fig. 4 and 5, respectively. 
Both cases show a good agreement of the model with the DNS data, indicating that 
the two models are equivalent in this case. It should be noted that the optimized C u 
values obtained for the u 2 u 2 model in the channel flow at lower Re (0.22), at high Re 
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FIGURE 4. y-distribution of the SGS-Reynolds shear stress in turbulent mixing 
layer. (Model with E = Eg/ /(AU ) 2 •) 



Figure 5. y-distributi on of the SGS-Reynolds shear stress in turbulent mixing 
layer, (model with E = u^u'^/iAU) 2 ) 

(0.23), and in the mixing layer (0.26) are very close to each other. This implies that 
the model constant of the SGS normal stress model is rather universal independent 
of the type of flow, whereas the optimized C v values for the E = Eg model were 
0.11, 0.10, 0.20, respectively. This is further indication of the potential strength of 
the normal-stress model. It is interesting to note that, particularly in the outer edge 
region of the mixing layer (y ~ 10 or ~ -10), the magnitude of normal component 
was the largest among the three components of the SGS turbulence fluctuations. A 
possible relationship of this phenomenon with the significant intermittency in these 
regions will be investigated in future work. 
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FIGURE 6. y-distribution of the SGS-Reynolds shear stress using Generalized nor- 
mal stress model (Channel at Rer = 790) 



FIGURE 7. y-distribution of the SGS-Reynolds shear stress using Generalized nor- 
mal stress model (Mixing layer) 

2.2 Generalized SGS normal- stress model 

Although it is shown t hat th e SGS-EVM based on the ‘major-gradient stress’ 
model (Eq. (4) with E — u^u^) shows a high correlation with the DNS data, this 
model does not preserve the tensoral invariance of the SGS Reynolds stresses. This 
drawback can be circumvented by generalizing this model as follows (Durbin 1991 
and Horiuti 1993): 


u\u) = S tj {\EG + \P) 


'tn 


du } 

dxi 


tjt 


dUi 

dxi 


( 7 ) 
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where 


= C t 


*V/2 

&G 


Ku'j, 


P = Ve„ 


du m 
1 dx, - 


This model is not identical to the conventional eddy viscosity model. In this model, 
the velocity scale is chosen as the component normal to the principal shear plane, 
for example, for u\u' 2 , the most dominant term in the RHS becomes nu tJ7 du i/x 2 . 

The generalized normal stress model is tested in the same manner as in the last 
section. The results for the high-J?e channel flow (C„ = 0.23) and the mixing 
layer ( C„ = 0.26) are shown in Figs. 6 and 7, respectively. It is evident that the 
generalized normal stress model shows a high correlation with the previous normal 
stress model. 

It should be pointed out that the generalization of the normal stress model with- 
out losing the tensorial invariance is not unique, e.g. terms such as 


—v. 


dui 


v, 


du ii 


( 8 ) 


can be added to Eq. (7). When these terms are in cluded , however, the term 
v x du\jdxi causes a large peak in the model for the u[ u 2 profile in the channel 
flow. For this reason, the terms in Eq. (8) were excluded when generalizing the 
normal stress model. 

2.3 Approximation method of the SGS turbulent energy 

To effectively use the model advocated in the previous sections, a model for the 
normal stresses is needed. We can either carry equations for the normal stresses 
or estimate the energy in the subgrid-scales from the energy in the large scales. 
In testing the scale-similarity model of Bardina (1983), Horiuti (1993) found good 
correlation between the model and the data. The model reads, 


Eg = Ck(ui - ui){u( - u t ), (9) 

u' 2 u' 2 = C N (u 2 - U 2 )(U 2 - W 2 ), 

where a constant different from unity was needed. It was pointed out that the 
optimized model constants Ck and Cn were not equal to unity because the scale- 
similarity model provides a partial estimate of the whole SGS fluctuations which 
resides in the vicinity of the cutoff- wave number (= tt/ A). The poor performance 
of the model when these coefficients are set equal to unity can be evidenced by 
the fact that in this case the SGS flatness parameter A becomes identically zero 
(purely two-dimensional state). We have optimized Ck/Cn for the low- J?e, the 
high-iZe channel, and the mixing layer flows and found 7.0/12.0, 7.0/ 9.0, and 9.0/ 
12.0 to be representative values for these flows. We note that they are slightly (but 
tolerably) sensitive to the type of flow fiel, and that they are generally close to each 
other. A representative comparison of the model prediction with the DNS data for 
Ujti i is shown for the high- Re channel flow and the mixing layer in Figs. 8 and 9, 
respectively. We find a good agreement with the DNS data for both flows. 
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FIGURE 8 . y-distribution of u\u\/u 2 r using the Bardina model (Channel at Re T — 
790) 



FIGURE 9. y-distribution of u\u\/{AU ) 2 using the Bardina model (Mixing layer) 

3. Future plans 

An ultimate goal of the present study is to develop an SGS model which yields 
good predictions of turbulent flows in a complex geometry. Of particular interest is 
the flow over a backward-facing step. In this case, while the flow is bounded by the 
walls, the internal mixing layer present in this flow plays a major role in setting the 
turbulence levels. In this work, a proper velocity scale for the SGS-EVM viscosity 
was determined for the fully developed channel and the mixing layer flows. In 
the channel, a clear advantage over more conventional treatments was shown by 
using the normal stress. It was also shown that the SGS normal stress is equally 
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useful as the total SGS turbulent energy for modeling in mixing layers. The model 
constant in the normal-stress model was found to be fairly independent of the type 
of flow. It was demonstrated that a generalized normal stress model can be used 
as an alternative method for the normal-stress model, and the tensorial invariance 
which is violated in the normal-stress model can be recovered. This generalized 
normal-stress model will be tested in a backward-facing step flow in both ‘a priori’ 
and 4 a posteriori’ manner in the future. Although the Bardina model constants Ck 
and Cn in Eq. (9) are rather consistent in three different flow fields, some variance 
was noticed. An attempt to determine these coefficients more accurately using the 
Dynamic scale model approach is currently underway. 
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